
#ifndef __libmesh_project_soln_h__
#define __libmesh_project_soln_h__

#include "projection_errors.h"

// libMesh includes
#include "mesh_function.h"

/**
 * This is the exact solution functor.
 */
class LibMeshProjection : public ProjectionBase
{
public:
  LibMeshProjection(bool verbose=false) ;

  virtual ~LibMeshProjection() ;

  virtual int apply(std::string output_file="") ;

  virtual std::string type() ;

//  virtual void init() ;

  class FunctionValue : public FunctionBase<Real>
  {
    public:
      FunctionValue(const LibMeshProjection& parent_obj, std::string p_sysname)
          : parent(&parent_obj), current_sys_name(p_sysname)
      {
        System &curr_sys = parent_obj.source->get_system(current_sys_name);

        unknown_names.resize(curr_sys.n_vars());
        for (unsigned int j = 0; j != curr_sys.n_vars(); ++j)
          unknown_names[j] = curr_sys.variable_name(j) ;
      }

      virtual ~FunctionValue() {}

      virtual Real operator() (const Point &p, const Real time=0.) ;

      virtual void operator() (const Point &p, const Real time, DenseVector<Real> &output) ;

      AutoPtr< FunctionBase<Real> > clone	()	 const ;

    private:
      const LibMeshProjection* parent ;
      std::vector<std::string> unknown_names;
      const std::string current_sys_name;
  };


/*
  Number fptr(const Point& p,
            const Parameters&,
            const std::string& libmesh_dbg_var(sys_name),
            const std::string& unknown_name) ;

  Gradient gptr(const Point& p,
              const Parameters&,
              const std::string& libmesh_dbg_var(sys_name),
              const std::string& unknown_name) ;
*/

private:

  mutable std::map<std::string, MeshFunction *> mesh_functions;

};

inline
LibMeshProjection::LibMeshProjection(bool is_verbose) : ProjectionBase(is_verbose)
{
}

inline
LibMeshProjection::~LibMeshProjection()
{
}


inline
std::string LibMeshProjection::type()
{
  return "LibMesh_Projection" ;
}


inline
Real LibMeshProjection::FunctionValue::operator() (const Point &p, const Real )
{
  libmesh_assert(parent->mesh_functions.count(unknown_names[0]));
  libmesh_assert(parent->mesh_functions[unknown_names[0]]);

  parent->projection_log->push("compute_value()");
  MeshFunction& meshfunc = *parent->mesh_functions[unknown_names[0]];
  Real value = meshfunc(p) ;
  parent->projection_log->pop("compute_value()");

  return value;
}


inline
void LibMeshProjection::FunctionValue::operator() (const Point &p, const Real , DenseVector<Real> &output)
{
  output.resize(unknown_names.size());
  for (unsigned int i=0; i < output.size(); ++i)
  {
    parent->projection_log->push("compute_value()");
    MeshFunction &meshfunc = *parent->mesh_functions[unknown_names[i]];
    output(i)= meshfunc(p) ;
    parent->projection_log->pop("compute_value()");
  }
}

inline
AutoPtr< FunctionBase<Real> > LibMeshProjection::FunctionValue::clone	()	 const
{
  return AutoPtr<FunctionBase<Real> > (new FunctionValue(*parent, current_sys_name)) ;
}

#endif // __libmesh_project_soln_h__
